数据包加载
knitr::opts_chunk$set(echo = TRUE)
# 检测和安装依赖包
package_list <- c("randomForest","ggplot2","pheatmap")
# 判断R包加载是否成功来决定是否安装后再加载
for(p in package_list){
if(!suppressWarnings(suppressMessages(require(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))){
install.packages(p)
suppressWarnings(suppressMessages(library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))
}
}
library(pheatmap)文件输入+自定义分组
# 1. 读取OTU表
otutab = read.table("C:/16S/result/randomforestclassification/otutab.txt", header=T, row.names= 1, sep="\t", comment.char = "", stringsAsFactors = F)
# 2. 读取物种注释
tax = read.table("C:/16S/result/randomforestclassification/taxonomy.txt", header=T, row.names= 1, sep="\t",comment.char = "", stringsAsFactors = F)
# 数据按实验设计手动筛选
metadata_all = read.table("C:/16S/result/randomforestclassification/metadata.txt", header=T, row.names=1, sep="\t")
SA = ""
SB = ""
SC = ""
Site_A = "S"
Site_B = "E"
Site_C = "G"
State_A = "H"
State_B = "D"
State_C = "I"
Sday_A = "pre"
Sday_B = "med"
Sday_C = "pre"
filetial1 = "rf.txt"
filetial2 = "rf.pdf"
sub_metadatanames = "sub_metadata"
rfcv_names = "rfcv"
importance_names = "importance"
top_f = "top_feautres"
prediction_name = "prediction_binary"
test_names = "test"
class_names = "new"
result_name = "result"
output_site = "C:/16S/result/randomforestclassification/"
result_file_names <- paste(result_name,Site_A,Site_B,"new",sep ="_")
result_file_names <- paste(output_site,result_file_names,sep ="")
sub_metadata_filename <- paste(sub_metadatanames,Site_A,Site_B,filetial1,sep ="_")
sub_metadata_filename <- paste(output_site,sub_metadata_filename,sep ="")
rfcv_filename_txt <- paste(class_names,rfcv_names,Site_A,Site_B,filetial1,sep ="_")
rfcv_filename_txt <- paste(output_site,rfcv_filename_txt,sep ="")
rfcv_filename_pdf <- paste(class_names,rfcv_names,Site_A,Site_B,filetial2,sep ="_")
rfcv_filename_pdf <- paste(output_site,rfcv_filename_pdf,sep ="")
importance_filename_txt <- paste(class_names,importance_names,Site_A,Site_B,filetial1,sep ="_")
importance_filename_txt <- paste(output_site,importance_filename_txt,sep ="")
top_feautre_filename <- paste(class_names,top_f,Site_A,Site_B,filetial2,sep ="_")
top_feautre_filename <- paste(output_site,top_feautre_filename,sep ="")
prediction_filename <- paste(class_names,prediction_name,Site_A,Site_B,filetial1,sep ="_")
prediction_filename <- paste(output_site,prediction_filename,sep ="")
testfilename1 <- paste(class_names,test_names,Site_A,filetial2,sep ="_")
testfilename1 <- paste(output_site,testfilename1,sep ="")
testfilename2 <- paste(class_names,test_names,Site_B,filetial2,sep ="_")
testfilename2 <- paste(output_site,testfilename2,sep ="")
print(result_file_names)## [1] "C:/16S/result/randomforestclassification/result_S_E_new"print(sub_metadata_filename)## [1] "C:/16S/result/randomforestclassification/sub_metadata_S_E_rf.txt"print(rfcv_filename_txt)## [1] "C:/16S/result/randomforestclassification/new_rfcv_S_E_rf.txt"print(rfcv_filename_pdf)## [1] "C:/16S/result/randomforestclassification/new_rfcv_S_E_rf.pdf"print(importance_filename_txt)## [1] "C:/16S/result/randomforestclassification/new_importance_S_E_rf.txt"print(top_feautre_filename)## [1] "C:/16S/result/randomforestclassification/new_top_feautres_S_E_rf.pdf"print(prediction_filename)## [1] "C:/16S/result/randomforestclassification/new_prediction_binary_S_E_rf.txt"print(testfilename1)## [1] "C:/16S/result/randomforestclassification/new_test_S_rf.pdf"print(testfilename2)## [1] "C:/16S/result/randomforestclassification/new_test_E_rf.pdf"# 筛选自己想要的分组以及内容
# sub_metadata1 = subset(metadata_all, State %in% c(State_A))
# sub_metadata1
sub_metadata2 = subset(metadata_all, Site_1 %in% c(Site_A,Site_B))
sub_metadata2# 实验设计与输入文件交叉筛选。
idx = rownames(sub_metadata2) %in% colnames(otutab)
sub_metadata3 = sub_metadata2[idx,]
sub_metadata3sub_otutab_tax = otutab[,rownames(sub_metadata3)]
sub_otutab_taxwrite.table(sub_metadata3,file = sub_metadata_filename,quote = F,sep = '\t', row.names = T, col.names = T)
# OTU丰度筛选阈值，默认0.1%，0为来筛选
thre = 0.02
# 输出文件名前缀
prefix = "tax_S_E_new"
# 生成各分类级汇总特征表
suppressWarnings(suppressMessages(library(amplicon)))
format2stamp(sub_otutab_tax, tax, thre, prefix)## Warning: 程辑包'dplyr'是用R版本4.3.2 来建造的now_dir <- getwd() # 当前工作目录
aim_dir <- "C:/16S/result/randomforestclassification"
# 创建目标文件夹（如果不存在）
if (!dir.exists(aim_dir)) {
dir.create(aim_dir, recursive = TRUE)
}
file_patterns <- c("tax_S_E_new_1Kingdom", "tax_S_E_new_2Phylum", "tax_S_E_new_3Class", "tax_S_E_new_4Order",
"tax_S_E_new_5Family", "tax_S_E_new_6Genus", "tax_S_E_new_7Species", "tax_S_E_new_8OTU0.02")
# 遍历文件模式并移动文件
sapply(file_patterns, function(pattern) {
files_to_move <- list.files(now_dir, pattern = pattern, full.names = TRUE)
sapply(files_to_move, function(file) {
file.copy(file, aim_dir, overwrite = TRUE)
file.remove(file)
})
})## tax_S_E_new_1Kingdom.C:/16S/All_figures_pipeline/tax_S_E_new_1Kingdom.txt
## TRUE
## tax_S_E_new_2Phylum.C:/16S/All_figures_pipeline/tax_S_E_new_2Phylum.txt
## TRUE
## tax_S_E_new_3Class.C:/16S/All_figures_pipeline/tax_S_E_new_3Class.txt
## TRUE
## tax_S_E_new_4Order.C:/16S/All_figures_pipeline/tax_S_E_new_4Order.txt
## TRUE
## tax_S_E_new_5Family.C:/16S/All_figures_pipeline/tax_S_E_new_5Family.txt
## TRUE
## tax_S_E_new_6Genus.C:/16S/All_figures_pipeline/tax_S_E_new_6Genus.txt
## TRUE
## tax_S_E_new_7Species.C:/16S/All_figures_pipeline/tax_S_E_new_7Species.txt
## TRUE
## tax_S_E_new_8OTU0.02.C:/16S/All_figures_pipeline/tax_S_E_new_8OTU0.02.txt
## TRUE# 打印成功信息
print("Files generated and moved to the target directory successfully.")## [1] "Files generated and moved to the target directory successfully."# 在指定目录生成tax_1-8共7个级别+OTU过滤文件分类级选择
先使用format2stamp.Rmd基于OTU表(otutab.txt)、物种注释(taxonomy.txt)和元数据(metadata.txt)筛选样本、高丰度特征，并分类汇总各分类级(tax_1-8)。 然后对各分类级进行准确性评估
# 读取实验设计、和物种分类文件
# 实验用的数据重命名
metadata_f = sub_metadata3
metadata_f# 设置训练集
train = subset(metadata_f, Site_1 %in% c("S","E"))
trainidx3 = sample(1:386, size = 200, replace = F)
# 注意sample后面的比例要根据上面文件自己调。size就是比例
# 选择的样本标记为TRUE，未选择的为FALSE
idx3 = 1:386 %in% idx3
#train2 = metadata_f[idx3,]
train2 = train[idx3,] # train2是根据实验设计在train里面又筛选的作为训练集的数据。测试集是剩余部分。
train2train2$Site_1 = as.factor(train2$Site_1)
summary(train2)## Filename Group_day Group_day1 Group_dayE
## Length:200 Length:200 Length:200 Length:200
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## Group_dayS Group StateSite StateSite_1
## Length:200 Length:200 Length:200 Length:200
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## State Site Site_1 Day
## Length:200 Length:200 E: 73 Min. : 0.000
## Class :character Class :character S:127 1st Qu.: 0.000
## Mode :character Mode :character Median : 3.000
## Mean : 3.435
## 3rd Qu.: 7.000
## Max. :10.000
## Day_num Day1 DayE Individual
## Min. : 0.000 Length:200 Length:200 Length:200
## 1st Qu.: 0.000 Class :character Class :character Class :character
## Median : 3.000 Mode :character Mode :character Mode :character
## Mean : 3.435
## 3rd Qu.: 7.000
## Max. :10.000# 筛选一个合适的分类等级
library(randomForest)
# 自己运行
for(i in c("2Phylum","3Class","4Order","5Family","6Genus","8OTU0.02")){
# i="5Family"
set.seed(0)
table_f = read.table(paste0("C:/16S/result/randomforestclassification/","tax_S_E_new_",i,".txt"),header = T, row.names = 1)
table_f = table_f[,rownames(train2)]
train2$Site_1 = as.factor(train2$Site_1)
rf = randomForest(t(table_f), train2$Site_1, importance=T, proximity=T, ntree = 5000)
print(i)
print(rf)
}## [1] "2Phylum"
##
## Call:
## randomForest(x = t(table_f), y = train2$Site_1, ntree = 5000, importance = T, proximity = T)
## Type of random forest: classification
## Number of trees: 5000
## No. of variables tried at each split: 3
##
## OOB estimate of error rate: 30%
## Confusion matrix:
## E S class.error
## E 29 44 0.6027397
## S 16 111 0.1259843
## [1] "3Class"
##
## Call:
## randomForest(x = t(table_f), y = train2$Site_1, ntree = 5000, importance = T, proximity = T)
## Type of random forest: classification
## Number of trees: 5000
## No. of variables tried at each split: 4
##
## OOB estimate of error rate: 29.5%
## Confusion matrix:
## E S class.error
## E 30 43 0.5890411
## S 16 111 0.1259843
## [1] "4Order"
##
## Call:
## randomForest(x = t(table_f), y = train2$Site_1, ntree = 5000, importance = T, proximity = T)
## Type of random forest: classification
## Number of trees: 5000
## No. of variables tried at each split: 5
##
## OOB estimate of error rate: 22.5%
## Confusion matrix:
## E S class.error
## E 34 39 0.53424658
## S 6 121 0.04724409
## [1] "5Family"
##
## Call:
## randomForest(x = t(table_f), y = train2$Site_1, ntree = 5000, importance = T, proximity = T)
## Type of random forest: classification
## Number of trees: 5000
## No. of variables tried at each split: 7
##
## OOB estimate of error rate: 23.5%
## Confusion matrix:
## E S class.error
## E 33 40 0.54794521
## S 7 120 0.05511811
## [1] "6Genus"
##
## Call:
## randomForest(x = t(table_f), y = train2$Site_1, ntree = 5000, importance = T, proximity = T)
## Type of random forest: classification
## Number of trees: 5000
## No. of variables tried at each split: 9
##
## OOB estimate of error rate: 18%
## Confusion matrix:
## E S class.error
## E 39 34 0.46575342
## S 2 125 0.01574803
## [1] "8OTU0.02"
##
## Call:
## randomForest(x = t(table_f), y = train2$Site_1, ntree = 5000, importance = T, proximity = T)
## Type of random forest: classification
## Number of trees: 5000
## No. of variables tried at each split: 23
##
## OOB estimate of error rate: 13.5%
## Confusion matrix:
## E S class.error
## E 47 26 0.356164384
## S 1 126 0.007874016# OOB estimate of error rate: 越小越好。最佳水平数据读取和统计
读取实验设计、Feature表，并进行数据筛选和交叉筛选
# 读取实验设计和物种分类文件
# 之前有过下面一样的步骤，但确保正确在重新运行一次。
train = subset(metadata_f, Site_1 %in% c("S","E"))
trainidx3 = sample(1:386, size = 200, replace = F)
# 注size=67，66比较小
# 选择的样本标记为TRUE，未选择的为FALSE
idx3 = 1:386 %in% idx3
#train2 = metadata_f[idx3,]
train2 = train[idx3,] # train2是根据实验设计在train里面又筛选的作为训练集的数据。测试集是剩余部分。
train2train2$Site_1 = as.factor(train2$Site_1)
summary(train2)## Filename Group_day Group_day1 Group_dayE
## Length:200 Length:200 Length:200 Length:200
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## Group_dayS Group StateSite StateSite_1
## Length:200 Length:200 Length:200 Length:200
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## State Site Site_1 Day Day_num
## Length:200 Length:200 E: 79 Min. : 0.00 Min. : 0.00
## Class :character Class :character S:121 1st Qu.: 0.00 1st Qu.: 0.00
## Mode :character Mode :character Median : 3.00 Median : 3.00
## Mean : 3.32 Mean : 3.32
## 3rd Qu.: 7.00 3rd Qu.: 7.00
## Max. :10.00 Max. :10.00
## Day1 DayE Individual
## Length:200 Length:200 Length:200
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
## # 根据之前筛选情况读取物种分类水平数据。读取对应水平特征表
table =read.table(paste0("C:/16S/result/randomforestclassification/","tax_S_E_new_6Genus.txt"),header = T, row.names = 1)
# train2重命名，方便后面区分。
metadata_train1 = train2
summary(metadata_train1)## Filename Group_day Group_day1 Group_dayE
## Length:200 Length:200 Length:200 Length:200
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## Group_dayS Group StateSite StateSite_1
## Length:200 Length:200 Length:200 Length:200
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## State Site Site_1 Day Day_num
## Length:200 Length:200 E: 79 Min. : 0.00 Min. : 0.00
## Class :character Class :character S:121 1st Qu.: 0.00 1st Qu.: 0.00
## Mode :character Mode :character Median : 3.00 Median : 3.00
## Mean : 3.32 Mean : 3.32
## 3rd Qu.: 7.00 3rd Qu.: 7.00
## Max. :10.00 Max. :10.00
## Day1 DayE Individual
## Length:200 Length:200 Length:200
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
## # 筛选OTU
idx = rownames(metadata_train1) %in% colnames(table)
metadata_train1 = metadata_train1[idx,]
otu_sub1 = table[, rownames(metadata_train1)]
otu_sub1dim(otu_sub1)## [1] 92 200选择最佳随机数(可选)
library(randomForest)
set.seed(77)
for (i in 0:9){
set.seed(i)
rf2 = randomForest(t(otu_sub1), metadata_train1$Site_1, importance=TRUE, proximity=TRUE, ntree = 5000)
print(i)
print(rf2)
}## [1] 0
##
## Call:
## randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000, importance = TRUE, proximity = TRUE)
## Type of random forest: classification
## Number of trees: 5000
## No. of variables tried at each split: 9
##
## OOB estimate of error rate: 17.5%
## Confusion matrix:
## E S class.error
## E 50 29 0.36708861
## S 6 115 0.04958678
## [1] 1
##
## Call:
## randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000, importance = TRUE, proximity = TRUE)
## Type of random forest: classification
## Number of trees: 5000
## No. of variables tried at each split: 9
##
## OOB estimate of error rate: 17%
## Confusion matrix:
## E S class.error
## E 52 27 0.34177215
## S 7 114 0.05785124
## [1] 2
##
## Call:
## randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000, importance = TRUE, proximity = TRUE)
## Type of random forest: classification
## Number of trees: 5000
## No. of variables tried at each split: 9
##
## OOB estimate of error rate: 18%
## Confusion matrix:
## E S class.error
## E 52 27 0.34177215
## S 9 112 0.07438017
## [1] 3
##
## Call:
## randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000, importance = TRUE, proximity = TRUE)
## Type of random forest: classification
## Number of trees: 5000
## No. of variables tried at each split: 9
##
## OOB estimate of error rate: 16.5%
## Confusion matrix:
## E S class.error
## E 52 27 0.34177215
## S 6 115 0.04958678
## [1] 4
##
## Call:
## randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000, importance = TRUE, proximity = TRUE)
## Type of random forest: classification
## Number of trees: 5000
## No. of variables tried at each split: 9
##
## OOB estimate of error rate: 17%
## Confusion matrix:
## E S class.error
## E 53 26 0.3291139
## S 8 113 0.0661157
## [1] 5
##
## Call:
## randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000, importance = TRUE, proximity = TRUE)
## Type of random forest: classification
## Number of trees: 5000
## No. of variables tried at each split: 9
##
## OOB estimate of error rate: 17.5%
## Confusion matrix:
## E S class.error
## E 51 28 0.35443038
## S 7 114 0.05785124
## [1] 6
##
## Call:
## randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000, importance = TRUE, proximity = TRUE)
## Type of random forest: classification
## Number of trees: 5000
## No. of variables tried at each split: 9
##
## OOB estimate of error rate: 18%
## Confusion matrix:
## E S class.error
## E 52 27 0.34177215
## S 9 112 0.07438017
## [1] 7
##
## Call:
## randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000, importance = TRUE, proximity = TRUE)
## Type of random forest: classification
## Number of trees: 5000
## No. of variables tried at each split: 9
##
## OOB estimate of error rate: 17.5%
## Confusion matrix:
## E S class.error
## E 53 26 0.32911392
## S 9 112 0.07438017
## [1] 8
##
## Call:
## randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000, importance = TRUE, proximity = TRUE)
## Type of random forest: classification
## Number of trees: 5000
## No. of variables tried at each split: 9
##
## OOB estimate of error rate: 16.5%
## Confusion matrix:
## E S class.error
## E 52 27 0.34177215
## S 6 115 0.04958678
## [1] 9
##
## Call:
## randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000, importance = TRUE, proximity = TRUE)
## Type of random forest: classification
## Number of trees: 5000
## No. of variables tried at each split: 9
##
## OOB estimate of error rate: 17%
## Confusion matrix:
## E S class.error
## E 53 26 0.3291139
## S 8 113 0.0661157随机森林分类
在确定的分类层级和最佳随机数下建模
library(randomForest)
set.seed(8)
rf3 = randomForest(t(otu_sub1), metadata_train1$Site_1, importance=TRUE, proximity=TRUE, ntree = 5000)
print(rf3)##
## Call:
## randomForest(x = t(otu_sub1), y = metadata_train1$Site_1, ntree = 5000, importance = TRUE, proximity = TRUE)
## Type of random forest: classification
## Number of trees: 5000
## No. of variables tried at each split: 9
##
## OOB estimate of error rate: 16.5%
## Confusion matrix:
## E S class.error
## E 52 27 0.34177215
## S 6 115 0.04958678交叉验证选择重要特征
set.seed(1) # 随机数据保证结果可重复，必须
# rfcv是随机森林交叉验证函数：Random Forest Cross Validation
result = rfcv(t(otu_sub1), metadata_train1$Site_1, cv.fold=5)
# 查看错误率表，31时错误率最低，为最佳模型
result$error.cv## 92 46 23 12 6 3 1
## 0.190 0.200 0.240 0.225 0.280 0.345 0.425# 绘制验证结果
with(result, plot(n.var, error.cv, log="x", type="o", lwd=2))
# 多次绘制
## 建立数据框保存多次结果
error.cv0 = data.frame(num = result$n.var, error.1 = result$error.cv)
## 指定随机数循环5次
for (i in 1:(1+9)){
print(i)
set.seed(i)
result= rfcv(t(otu_sub1), metadata_train1$Site_1, cv.fold=5) # scale = "log", step = 0.9
error.cv0 = cbind(error.cv0, result$error.cv)
}## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10error.cv0 绘制交叉验证曲线
# 提取x轴标签
n.var = error.cv0$num
# 提取y轴数据+标签
error.cv = error.cv0[,2:6]
colnames(error.cv) = paste('err',1:5,sep='.')
# 添加均值
err.mean = apply(error.cv,1,mean)
# 合并新的数据库，x+error+mean
allerr = data.frame(num=n.var,err.mean=err.mean,error.cv)
# number of otus selected 人为在图中观察的结果，30几乎为最低，且数量可接受
optimal = 20
# 图1：机器学习结果交叉验证图，选择Top features
# 图中 + 5条灰色拆线+1条黑色均值拆线+一条最优垂线+X轴对数变换
write.table(allerr, file = rfcv_filename_txt, sep = "\t", quote = F, row.names = T, col.names = T)
p = ggplot() + # 开始绘图
geom_line(aes(x = allerr$num, y = allerr$err.1), colour = 'grey') + # 5次验证灰线
geom_line(aes(x = allerr$num, y = allerr$err.2), colour = 'grey') +
geom_line(aes(x = allerr$num, y = allerr$err.3), colour = 'grey') +
geom_line(aes(x = allerr$num, y = allerr$err.4), colour = 'grey') +
geom_line(aes(x = allerr$num, y = allerr$err.5), colour = 'grey') +
geom_line(aes(x = allerr$num, y = allerr$err.mean), colour = 'black') + # 均值黑线
geom_vline(xintercept = optimal, colour='black', lwd=0.36, linetype="dashed") + # 最优垂线
coord_trans(x = "log2") + # X轴对数变换和刻度
scale_x_continuous(breaks = c(1, 2, 5, 10, 20, 30, 50, 100, 200)) + # , max(allerr$num)
labs(title=paste('Training set (n = ', dim(t(otu_sub1))[1],')', sep = ''),
x='Number of families ', y='Cross-validation error rate') +
annotate("text", x = optimal, y = max(allerr$err.mean), label=paste("optimal = ", optimal, sep="")) + theme_bw()
p 
ggsave(p, file = rfcv_filename_pdf, width = 180, height = 120, unit = 'mm')特征重要性可视化
## 预览和保存特征贡献度
imp= as.data.frame(rf3$importance)
imp = imp[order(imp$MeanDecreaseAccuracy, decreasing = T),]
# bbb = head(imp,n=optimal)
write.table(imp,file = importance_filename_txt,quote = F,sep = '\t', row.names = T, col.names = T)
# 读取所有feature贡献度
imp2 = read.table(importance_filename_txt, header=T, row.names= 1, sep="\t")
# 分析选择top20分组效果最好，参数显示数量
imp2 = head(imp2, n = optimal)
imp2 = imp2[order(imp2$MeanDecreaseAccuracy, decreasing = F),]
# 简化全名，去掉界
imp2$Family = gsub("Bacteria\\|","",rownames(imp2))
# 添加门用于着色(删除竖线后面全部)
imp2$Phylum = gsub("\\|.*","",imp2$Family)
# 设置顺序
imp2$Family = factor(imp2$Family, levels = imp2$Family)
# 图2. 绘制物种类型种重要性柱状图
p = ggplot(imp2, aes(x = Family, y = MeanDecreaseAccuracy, fill = Phylum)) +
geom_bar(stat = "identity") +
coord_flip() + theme_bw()
p
ggsave(top_feautre_filename, p, width=200*2.5, height=59*2, unit='mm')
# 名称不一定唯一，需要手动修改测试集独立验证
如果第一地点数据量足够大，可以取出1/2到1/3进行同一地点的独立验证。方法相同。 筛选测序集样品
# 测试集选择
test = train[!idx3,]
test#test = metadata_f[!idx3,]
metadata_test = test
summary(metadata_test)## Filename Group_day Group_day1 Group_dayE
## Length:186 Length:186 Length:186 Length:186
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## Group_dayS Group StateSite StateSite_1
## Length:186 Length:186 Length:186 Length:186
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## State Site Site_1 Day
## Length:186 Length:186 Length:186 Min. : 0.000
## Class :character Class :character Class :character 1st Qu.: 0.000
## Mode :character Mode :character Mode :character Median : 3.000
## Mean : 3.312
## 3rd Qu.: 7.000
## Max. :10.000
## Day_num Day1 DayE Individual
## Min. : 0.000 Length:186 Length:186 Length:186
## 1st Qu.: 0.000 Class :character Class :character Class :character
## Median : 3.000 Mode :character Mode :character Mode :character
## Mean : 3.312
## 3rd Qu.: 7.000
## Max. :10.000idx4 = rownames(metadata_test) %in% colnames(table)
metadata_test = metadata_test[idx4,]
otu_sub2 = table[,rownames(metadata_test)]
# 转置，并添加分组信息
otutab_t2 = as.data.frame(t(otu_sub2))
otutab_t2$Site_1 = metadata_test[rownames(otutab_t2),]$Site_1基于训练集随机森林模型验证
set.seed(315)
otutab.pred = predict(rf3, t(otu_sub2) )
pre_tab = table(observed=otutab_t2[,"Site_1"],
predicted=otutab.pred)
pre_tab## predicted
## observed E S
## E 54 22
## S 10 100可视化验证结果
# 整理样本原始分组和预测分类
predict = data.frame(group = otutab_t2[,"Site_1"], predicted=otutab.pred)
# 保存预测结果表
write.table("SampleID\t", file=prediction_filename,append = F, quote = F, eol = "", row.names = F, col.names = F)
write.table(predict, file = prediction_filename,append = T, quote = F, row.names = T, col.names = T, sep = "\t")## Warning in write.table(predict, file = prediction_filename, append = T, :
## 给文件加列名# 转换为数值可视化
# 预测准确标为1，错误标为0
predict$result = ifelse(predict$group == predict$predicted, 1, 0)
# IND=1, TEJ=2
predict$predict = ifelse(predict$predicted == "E", 1, 2)
# Set sample number in each row
column = 10
low_to_high_colors <- colorRampPalette(c("#B1DDE1", "#F8C7BD"))(100)
# Draw prediction result
AA1 = predict[predict$group=="S",]$predict
length(AA1)## [1] 110row = round(length(AA1)/column + 0.5)
row## [1] 12i = column * row - length(AA1)
AA1 = c(AA1, rep(NA, i))
matrix1 = matrix(AA1, nrow = row, ncol = column, byrow = T)
pheatmap(matrix1, cluster_rows = F, cluster_cols = F, cellwidth = 15, cellheight = 12, color = low_to_high_colors, border_color = "#EFEFEF")
pheatmap(matrix1, cluster_rows = F, cluster_cols = F, cellwidth = 18, cellheight = 18,filename = testfilename1, color = low_to_high_colors, border_color = "#EFEFEF")
# Draw another prediction result
BB1 = predict[predict$group=="E",]$predict
length(BB1)## [1] 76row = round(length(BB1)/column + 0.5)
i = column * row - length(BB1)
BB1 = c(BB1, rep(NA, i))
matrix2 = matrix(BB1, nrow = row, ncol = column, byrow = T)
pheatmap(matrix2, cluster_rows = F, cluster_cols = F, cellwidth = 15, cellheight = 12, color = low_to_high_colors, border_color = "#EFEFEF")
# 保存图片
pheatmap(matrix2, cluster_rows = F, cluster_cols = F, cellwidth = 18, cellheight = 18, filename = testfilename2, color = low_to_high_colors, border_color = "#EFEFEF")
## 将生成的文件放入指定文件夹
source_dir <- "C:/16S/result/randomforestclassification"
target_dir <- result_file_names
if (!dir.exists(target_dir)) {
dir.create(target_dir)
}
# 选取含特定字符的文件（例如后缀为 csv 的所有文件）
tobeCopy1 <- list.files(source_dir, pattern = "*_rf.txt", full.names = TRUE)
tobeCopy2 <- list.files(source_dir, pattern = "*_rf.pdf", full.names = TRUE)
# 复制选中文件到目标文件夹
sapply(tobeCopy1, function(x) { file.copy(x, target_dir, overwrite = TRUE) })## C:/16S/result/randomforestclassification/new_importance_S_E_rf.txt
## TRUE
## C:/16S/result/randomforestclassification/new_prediction_binary_S_E_rf.txt
## TRUE
## C:/16S/result/randomforestclassification/new_rfcv_S_E_rf.txt
## TRUE
## C:/16S/result/randomforestclassification/sub_metadata_S_E_rf.txt
## TRUEsapply(tobeCopy2, function(x) { file.copy(x, target_dir, overwrite = TRUE) })## C:/16S/result/randomforestclassification/new_rfcv_S_E_rf.pdf
## TRUE
## C:/16S/result/randomforestclassification/new_test_E_rf.pdf
## TRUE
## C:/16S/result/randomforestclassification/new_test_S_rf.pdf
## TRUE
## C:/16S/result/randomforestclassification/new_top_feautres_S_E_rf.pdf
## TRUE# 删除源文件夹中的已复制文件
file.remove(tobeCopy1)## [1] TRUE TRUE TRUE TRUEfile.remove(tobeCopy2)## [1] TRUE TRUE TRUE TRUE# 确认文件已经复制和删除
print("Files copied and deleted successfully.")## [1] "Files copied and deleted successfully."使用此脚本，请引用下文：
If used this script, please cited:
Yong-Xin Liu, Lei Chen, Tengfei Ma, Xiaofang Li, Maosheng Zheng, Xin Zhou, Liang Chen, Xubo Qian, Jiao Xi, Hongye Lu, Huiluo Cao, Xiaoya Ma, Bian Bian, Pengfan Zhang, Jiqiu Wu, Ren-You Gan, Baolei Jia, Linyang Sun, Zhicheng Ju, Yunyun Gao, Tao Wen, Tong Chen. 2023. EasyAmplicon: An easy-to-use, open-source, reproducible, and community-based pipeline for amplicon data analysis in microbiome research. iMeta 2: e83. https://doi.org/10.1002/imt2.83
Copyright 2016-2023 Yong-Xin Liu liuyongxin@caas.cn, Tao Wen taowen@njau.edu.cn, Tong Chen chent@nrc.ac.cn